Run brm of each transformed variable by “Years since infection”
First, set instructions for sampler
n_cores <- detectCores() # this will determine how many cores your computer has so that it can use all of its processing power
n_chains <- 2
n_iter <- 8000
n_warmup <- 2000
#-Also add a function to make it easier when we're working with logit-transformed data-#
logit <- function(mu) {
logit_mu <- log(mu/ (1 - mu))
return(logit_mu)
}
#-Set inits-#
init_func <- function(chain_id = 1) {
list ( Intercept = runif(1, .01, .1),
b = runif(1, .001, .01),
phi = runif(1, .001, .01),
sd = runif(1, .001, .01),
cor = runif(1, 0.001, .01),
zi = runif(1, 0.001, .01))
}
# Create different inits for each chain
init_list <- vector(mode = "list", length = n_chains)
for(i in 1:n_chains)
{
init_list[[i]] <- init_func(chain_id = i)
}
CCA model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_CalgCov)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

# subset just the Black Point site data because this will be reference level
BlackPointData <- tcrmp_meta_pctcov %>%
filter(Site=="Black Point")
#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_CalgCov) #-we'll work with the zero-inflated beta

#define model
cca_mod <- bf(pctCov_CalgCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
##### Get priors
get_prior(cca_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## (flat) b zi
## (flat) b yearsSinceInf zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd Site zi 0
## student_t(3, 0, 2.5) sd Intercept Site zi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site zi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_Calg <- BlackPointData %>%
filter(pctCov_CalgCov> 0)
logit(mean(BPD_Calg$pctCov_CalgCov)) #-5.2
## [1] -5.213258
##### Set priors
cca_prior <- c(prior(normal(-5.2, 1), class = "Intercept"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(normal(0, 1), class = "b"),
prior(normal(0, 1), class = "b", dpar = "zi"),
prior(lkj(2), class = "cor"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"),
prior(exponential(1), class = "sd", dpar = "zi"))
##### Uncomment below to run model
# cca_brm <- brm(cca_mod,
# data = tcrmp_meta_pctcov,
# prior = cca_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Save model and/or read in saved model
#saveRDS(cca_brm, "data/outputs/cca_brm.RDS")
cca_brm <- readRDS("data/outputs/cca_brm.RDS")
##### Look at model
cca_brm
## Family: zero_inflated_beta
## Links: mu = logit; phi = log; zi = logit
## Formula: pctCov_CalgCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.75 0.10 0.58 0.98 1.00
## sd(yearsSinceInf) 0.24 0.07 0.13 0.39 1.00
## sd(phi_Intercept) 0.92 0.15 0.66 1.26 1.00
## sd(phi_yearsSinceInf) 0.34 0.19 0.03 0.74 1.00
## sd(zi_Intercept) 2.32 0.40 1.67 3.25 1.00
## sd(zi_yearsSinceInf) 0.57 0.34 0.05 1.37 1.00
## cor(Intercept,yearsSinceInf) 0.37 0.23 -0.11 0.77 1.00
## cor(phi_Intercept,phi_yearsSinceInf) 0.22 0.34 -0.47 0.82 1.00
## cor(zi_Intercept,zi_yearsSinceInf) -0.18 0.36 -0.82 0.54 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2770 5141
## sd(yearsSinceInf) 5386 7242
## sd(phi_Intercept) 4028 7082
## sd(phi_yearsSinceInf) 2555 3744
## sd(zi_Intercept) 3561 6221
## sd(zi_yearsSinceInf) 2256 3399
## cor(Intercept,yearsSinceInf) 7530 8289
## cor(phi_Intercept,phi_yearsSinceInf) 9830 7529
## cor(zi_Intercept,zi_yearsSinceInf) 9834 7959
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -4.18 0.13 -4.44 -3.92 1.00 1655 2835
## phi_Intercept 4.92 0.18 4.56 5.27 1.00 3197 5215
## zi_Intercept -2.01 0.44 -2.92 -1.17 1.00 2244 3352
## yearsSinceInf -0.30 0.07 -0.43 -0.17 1.00 6367 7854
## phi_yearsSinceInf 0.38 0.14 0.11 0.65 1.00 8082 6806
## zi_yearsSinceInf 0.78 0.22 0.38 1.24 1.00 7904 7832
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(cca_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(cca_brm, type = "boxplot", ndraws = 6)

pp_check(cca_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(cca_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_CalgCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Cyanobacteria model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_CyanCov)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_CyanCov) #-we'll work with the zero-inflated beta

#define model
cyano_mod <- bf(pctCov_CyanCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
##### Get priors
get_prior(cyano_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## (flat) b zi
## (flat) b yearsSinceInf zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd Site zi 0
## student_t(3, 0, 2.5) sd Intercept Site zi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site zi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_Cyano <- BlackPointData %>%
filter(pctCov_CyanCov> 0)
logit(mean(BPD_Cyano$pctCov_CyanCov)) #-3.04
## [1] -3.043492
##### Set priors
cyano_prior <- c(prior(normal(-3, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(normal(0, 1), class = "b", dpar = "zi"),
prior(exponential(1), class = "sd", dpar = "zi"),
prior(lkj(2), class = "cor"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"))
##### Uncomment below to run model
# cyano_brm <- brm(cyano_mod,
# data = tcrmp_meta_pctcov,
# prior = cyano_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Save model and/or read in saved model
#saveRDS(cyano_brm, "data/outputs/cyanobacteria_brm.RDS")
cyano_brm <- readRDS("data/outputs/cyanobacteria_brm.RDS")
##### Look at model
cyano_brm
## Warning: There were 10 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Family: zero_inflated_beta
## Links: mu = logit; phi = log; zi = logit
## Formula: pctCov_CyanCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.68 0.10 0.51 0.91 1.00
## sd(yearsSinceInf) 0.71 0.13 0.48 0.99 1.00
## sd(phi_Intercept) 0.88 0.13 0.66 1.17 1.00
## sd(phi_yearsSinceInf) 0.47 0.27 0.03 1.02 1.00
## sd(zi_Intercept) 0.93 0.16 0.65 1.27 1.00
## sd(zi_yearsSinceInf) 1.63 0.49 0.84 2.72 1.00
## cor(Intercept,yearsSinceInf) 0.05 0.22 -0.37 0.47 1.00
## cor(phi_Intercept,phi_yearsSinceInf) -0.43 0.34 -0.91 0.40 1.00
## cor(zi_Intercept,zi_yearsSinceInf) 0.07 0.22 -0.35 0.49 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2414 5154
## sd(yearsSinceInf) 3001 6003
## sd(phi_Intercept) 2954 4270
## sd(phi_yearsSinceInf) 1070 2075
## sd(zi_Intercept) 4594 6293
## sd(zi_yearsSinceInf) 3536 5306
## cor(Intercept,yearsSinceInf) 2706 4261
## cor(phi_Intercept,phi_yearsSinceInf) 4957 5097
## cor(zi_Intercept,zi_yearsSinceInf) 3396 5732
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.57 0.12 -3.82 -3.32 1.00 1428 3003
## phi_Intercept 3.61 0.16 3.29 3.93 1.00 2005 3293
## zi_Intercept -1.16 0.18 -1.53 -0.80 1.00 3122 6107
## yearsSinceInf 0.15 0.14 -0.14 0.43 1.00 2242 4272
## phi_yearsSinceInf 0.19 0.16 -0.10 0.55 1.00 2703 4353
## zi_yearsSinceInf -0.86 0.38 -1.68 -0.17 1.00 4010 5783
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# posterior predictive checks
pp_check(cyano_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(cyano_brm, type = "boxplot", ndraws = 6)

pp_check(cyano_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(cyano_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_CyanCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Fire coral model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_HydrozoanCoral)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_HydrozoanCoral) #-we'll work with zero-inflated beta

#define model
fireCoral_mod <- bf(pctCov_HydrozoanCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
##### Get priors
get_prior(fireCoral_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## (flat) b zi
## (flat) b yearsSinceInf zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd Site zi 0
## student_t(3, 0, 2.5) sd Intercept Site zi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site zi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_fireCoral <- BlackPointData %>%
filter(pctCov_HydrozoanCoral> 0)
logit(mean(BPD_fireCoral$pctCov_HydrozoanCoral)) #-5.5
## [1] -5.450887
##### Set priors
fireCoral_prior <- c(prior(normal(-5.5, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(lkj(2), class = "cor"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(normal(0, 1), class = "b", dpar = "zi"),
prior(exponential(1), class = "sd", dpar = "zi"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"))
##### Uncomment to run model
# fireCoral_brm <- brm(fireCoral_mod,
# data = tcrmp_meta_pctcov,
# prior = fireCoral_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Uncomment below to save model and/or read in saved model
#saveRDS(fireCoral_brm, "data/outputs/fireCoral_brm.RDS")
fireCoral_brm <- readRDS("data/outputs/fireCoral_brm.RDS")
##### Look at model
fireCoral_brm
## Warning: There were 1 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Family: zero_inflated_beta
## Links: mu = logit; phi = log; zi = logit
## Formula: pctCov_HydrozoanCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.36 0.06 0.26 0.49 1.00
## sd(yearsSinceInf) 0.07 0.05 0.00 0.19 1.00
## sd(phi_Intercept) 0.79 0.14 0.56 1.11 1.00
## sd(phi_yearsSinceInf) 0.17 0.14 0.01 0.53 1.00
## sd(zi_Intercept) 1.55 0.24 1.15 2.08 1.00
## sd(zi_yearsSinceInf) 0.32 0.23 0.01 0.86 1.00
## cor(Intercept,yearsSinceInf) 0.07 0.42 -0.73 0.81 1.00
## cor(phi_Intercept,phi_yearsSinceInf) -0.05 0.40 -0.78 0.75 1.00
## cor(zi_Intercept,zi_yearsSinceInf) 0.05 0.41 -0.73 0.80 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 4345 7094
## sd(yearsSinceInf) 4756 5939
## sd(phi_Intercept) 4246 6768
## sd(phi_yearsSinceInf) 5087 6581
## sd(zi_Intercept) 3683 6402
## sd(zi_yearsSinceInf) 2769 5093
## cor(Intercept,yearsSinceInf) 13087 8056
## cor(phi_Intercept,phi_yearsSinceInf) 12654 8736
## cor(zi_Intercept,zi_yearsSinceInf) 13846 7685
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -5.21 0.08 -5.36 -5.06 1.00 3098 5274
## phi_Intercept 6.37 0.19 6.00 6.72 1.00 4121 5885
## zi_Intercept 0.98 0.29 0.41 1.55 1.00 2127 4044
## yearsSinceInf 0.08 0.05 -0.03 0.19 1.00 9591 9106
## phi_yearsSinceInf -0.10 0.15 -0.41 0.20 1.00 9744 8595
## zi_yearsSinceInf -0.22 0.15 -0.50 0.11 1.00 8079 5170
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# posterior predictive checks
pp_check(fireCoral_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(fireCoral_brm, type = "boxplot", ndraws = 6)
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?
## Notch went outside hinges
## ℹ Do you want `notch = FALSE`?

pp_check(fireCoral_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(fireCoral_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_HydrozoanCoral*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Gorgonian model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_GorgCov)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_GorgCov) #-we'll work with zero-inflated beta

#define model
gorgo_mod <- bf(pctCov_GorgCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
##### Get priors
get_prior(gorgo_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## (flat) b zi
## (flat) b yearsSinceInf zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd Site zi 0
## student_t(3, 0, 2.5) sd Intercept Site zi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site zi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_gorgo <- BlackPointData %>%
filter(pctCov_GorgCov> 0)
logit(mean(BPD_gorgo$pctCov_GorgCov)) #-3.6
## [1] -3.562578
##### Set priors
gorgo_prior <- c(prior(normal(-3.6, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(lkj(2), class = "cor"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(normal(0, 1), class = "b", dpar = "zi"),
prior(exponential(1), class = "sd", dpar = "zi"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"))
##### Uncomment to run model
# gorgo_brm <- brm(gorgo_mod,
# data = tcrmp_meta_pctcov,
# prior = gorgo_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Save model and/or read in saved model
#saveRDS(gorgo_brm, "data/outputs/gorgonian_brm.RDS")
gorgo_brm<- readRDS("data/outputs/gorgonian_brm.RDS")
##### Look at model
gorgo_brm
## Family: zero_inflated_beta
## Links: mu = logit; phi = log; zi = logit
## Formula: pctCov_GorgCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.99 0.13 0.78 1.27 1.00
## sd(yearsSinceInf) 0.09 0.05 0.01 0.21 1.00
## sd(phi_Intercept) 0.99 0.15 0.73 1.32 1.00
## sd(phi_yearsSinceInf) 0.14 0.10 0.01 0.39 1.00
## sd(zi_Intercept) 2.69 0.47 1.92 3.75 1.00
## sd(zi_yearsSinceInf) 0.39 0.25 0.02 0.94 1.00
## cor(Intercept,yearsSinceInf) -0.23 0.35 -0.83 0.52 1.00
## cor(phi_Intercept,phi_yearsSinceInf) 0.08 0.41 -0.72 0.81 1.00
## cor(zi_Intercept,zi_yearsSinceInf) -0.32 0.39 -0.90 0.58 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2856 5202
## sd(yearsSinceInf) 4069 5645
## sd(phi_Intercept) 4178 6217
## sd(phi_yearsSinceInf) 5453 7333
## sd(zi_Intercept) 3598 6354
## sd(zi_yearsSinceInf) 4396 5902
## cor(Intercept,yearsSinceInf) 16965 8601
## cor(phi_Intercept,phi_yearsSinceInf) 19649 8670
## cor(zi_Intercept,zi_yearsSinceInf) 14832 9383
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.85 0.18 -4.21 -3.51 1.00 1305 2642
## phi_Intercept 4.66 0.18 4.30 5.03 1.00 2973 5193
## zi_Intercept -2.52 0.52 -3.59 -1.54 1.00 2543 4046
## yearsSinceInf 0.01 0.04 -0.08 0.09 1.00 15958 9584
## phi_yearsSinceInf 0.03 0.10 -0.15 0.22 1.00 13837 9178
## zi_yearsSinceInf 0.08 0.22 -0.33 0.54 1.00 11903 8936
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(gorgo_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(gorgo_brm, type = "boxplot", ndraws = 6)

pp_check(gorgo_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(gorgo_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_GorgCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Macroalgae model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_MacaCov)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_MacaCov) #-we'll work with beta

#define model
macro_mod <- bf(pctCov_MacaCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = Beta(link = "logit", link_phi = "log"))
##### Get priors
get_prior(macro_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_macro <- BlackPointData %>%
filter(pctCov_MacaCov> 0)
logit(mean(BPD_macro$pctCov_MacaCov)) #-1
## [1] -0.9637347
##### Set priors
macro_prior <- c(prior(normal(-1, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(lkj(2), class = "cor"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"))
##### Uncomment to run model
# macro_brm <- brm(macro_mod,
# data = tcrmp_meta_pctcov,
# prior = macro_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Save model and/or read in saved model
#saveRDS(macro_brm, "data/outputs/macroalgae_brm.RDS")
macro_brm<- readRDS("data/outputs/macroalgae_brm.RDS")
##### Look at model
macro_brm
## Family: beta
## Links: mu = logit; phi = log
## Formula: pctCov_MacaCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.89 0.11 0.69 1.14 1.00
## sd(yearsSinceInf) 0.35 0.07 0.24 0.49 1.00
## sd(phi_Intercept) 0.56 0.09 0.41 0.75 1.00
## sd(phi_yearsSinceInf) 0.27 0.14 0.02 0.55 1.00
## cor(Intercept,yearsSinceInf) -0.61 0.13 -0.82 -0.30 1.00
## cor(phi_Intercept,phi_yearsSinceInf) -0.26 0.31 -0.77 0.47 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2026 3900
## sd(yearsSinceInf) 3924 6727
## sd(phi_Intercept) 4654 6889
## sd(phi_yearsSinceInf) 2417 3204
## cor(Intercept,yearsSinceInf) 4609 6972
## cor(phi_Intercept,phi_yearsSinceInf) 10667 7094
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.71 0.15 -1.02 -0.41 1.00 797 1435
## phi_Intercept 3.11 0.11 2.90 3.32 1.00 3699 5961
## yearsSinceInf 0.12 0.07 -0.01 0.27 1.00 1711 3440
## phi_yearsSinceInf 0.48 0.10 0.29 0.69 1.00 9179 8812
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(macro_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(macro_brm, type = "boxplot", ndraws = 6)

pp_check(macro_brm, type = "dens_overlay", ndraws = 100) #looks good

y <- posterior_predict(macro_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_MacaCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

SCTLD-resistant coral model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_ResistantCoral)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_ResistantCoral) #-we'll work with zero-inflated beta

#define model
resistantCoral_mod <- bf(pctCov_ResistantCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
##### Get priors
get_prior(resistantCoral_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## (flat) b zi
## (flat) b yearsSinceInf zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd Site zi 0
## student_t(3, 0, 2.5) sd Intercept Site zi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site zi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_resistantCoral <- BlackPointData %>%
filter(pctCov_ResistantCoral> 0)
logit(mean(BPD_resistantCoral$pctCov_ResistantCoral)) #-2.8
## [1] -2.807501
##### Set priors
resistantCoral_prior <- c(prior(normal(-2.8, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(lkj(2), class = "cor"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(normal(0, 1), class = "b", dpar = "zi"),
prior(exponential(1), class = "sd", dpar = "zi"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"))
##### Uncomment to run model
# resistantCoral_brm <- brm(resistantCoral_mod,
# data = tcrmp_meta_pctcov,
# prior = resistantCoral_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Uncomment below to save model and/or read in saved model
#saveRDS(resistantCoral_brm, "data/outputs/resistantCoral_brm.RDS")
resistantCoral_brm <- readRDS("data/outputs/resistantCoral_brm.RDS")
##### Look at model
resistantCoral_brm
## Family: zero_inflated_beta
## Links: mu = logit; phi = log; zi = logit
## Formula: pctCov_ResistantCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.71 0.09 0.56 0.91 1.00
## sd(yearsSinceInf) 0.05 0.04 0.00 0.16 1.00
## sd(phi_Intercept) 0.73 0.10 0.55 0.96 1.00
## sd(phi_yearsSinceInf) 0.11 0.08 0.00 0.31 1.00
## sd(zi_Intercept) 2.06 0.49 1.28 3.20 1.00
## sd(zi_yearsSinceInf) 0.31 0.28 0.01 1.04 1.00
## cor(Intercept,yearsSinceInf) -0.21 0.41 -0.84 0.69 1.00
## cor(phi_Intercept,phi_yearsSinceInf) -0.16 0.43 -0.87 0.71 1.00
## cor(zi_Intercept,zi_yearsSinceInf) 0.07 0.44 -0.78 0.84 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1870 4118
## sd(yearsSinceInf) 3458 5364
## sd(phi_Intercept) 3240 5804
## sd(phi_yearsSinceInf) 4628 5224
## sd(zi_Intercept) 3423 5595
## sd(zi_yearsSinceInf) 5332 5972
## cor(Intercept,yearsSinceInf) 8925 7147
## cor(phi_Intercept,phi_yearsSinceInf) 10718 8626
## cor(zi_Intercept,zi_yearsSinceInf) 9627 8870
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.59 0.12 -3.84 -3.35 1.00 903 1894
## phi_Intercept 4.69 0.14 4.42 4.95 1.00 2223 4119
## zi_Intercept -4.48 0.55 -5.69 -3.54 1.00 3987 5574
## yearsSinceInf 0.02 0.03 -0.05 0.08 1.00 8267 7681
## phi_yearsSinceInf 0.05 0.08 -0.11 0.21 1.00 10306 9134
## zi_yearsSinceInf -0.15 0.33 -0.86 0.45 1.00 8257 5633
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(resistantCoral_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(resistantCoral_brm, type = "boxplot", ndraws = 6)

pp_check(resistantCoral_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(resistantCoral_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_ResistantCoral*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

SCTLD-susceptible coral model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_SuscCoral)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_SuscCoral) #-we'll work with zero-inflated beta

#define model
SuscCoral_mod <- bf(pctCov_SuscCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
##### Get priors
get_prior(SuscCoral_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## (flat) b zi
## (flat) b yearsSinceInf zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd Site zi 0
## student_t(3, 0, 2.5) sd Intercept Site zi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site zi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_SuscCoral <- BlackPointData %>%
filter(pctCov_SuscCoral> 0)
logit(mean(BPD_SuscCoral$pctCov_SuscCoral)) #-2.4
## [1] -2.439212
##### Set priors
SuscCoral_prior <- c(prior(normal(-2.4, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(lkj(2), class = "cor"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(normal(0, 1), class = "b", dpar = "zi"),
prior(exponential(1), class = "sd", dpar = "zi"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"))
##### Uncomment to run model
# SuscCoral_brm <- brm(SuscCoral_mod,
# data = tcrmp_meta_pctcov,
# prior = SuscCoral_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Save model and/or read in saved model
#saveRDS(SuscCoral_brm, "data/outputs/susceptibleCoral_brm.RDS")
SuscCoral_brm <- readRDS("data/outputs/susceptibleCoral_brm.RDS")
##### Look at model
SuscCoral_brm
## Family: zero_inflated_beta
## Links: mu = logit; phi = log; zi = logit
## Formula: pctCov_SuscCoral ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 1.12 0.15 0.88 1.46 1.00
## sd(yearsSinceInf) 0.30 0.07 0.18 0.47 1.00
## sd(phi_Intercept) 0.73 0.11 0.54 0.96 1.00
## sd(phi_yearsSinceInf) 0.10 0.08 0.00 0.29 1.00
## sd(zi_Intercept) 3.08 0.65 2.03 4.56 1.00
## sd(zi_yearsSinceInf) 0.60 0.51 0.02 1.90 1.00
## cor(Intercept,yearsSinceInf) 0.05 0.23 -0.39 0.48 1.00
## cor(phi_Intercept,phi_yearsSinceInf) -0.05 0.43 -0.81 0.77 1.00
## cor(zi_Intercept,zi_yearsSinceInf) 0.23 0.43 -0.69 0.88 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1651 2979
## sd(yearsSinceInf) 3657 6269
## sd(phi_Intercept) 3920 6756
## sd(phi_yearsSinceInf) 5740 6482
## sd(zi_Intercept) 4207 6759
## sd(zi_yearsSinceInf) 3288 6503
## cor(Intercept,yearsSinceInf) 5678 8110
## cor(phi_Intercept,phi_yearsSinceInf) 15435 8051
## cor(zi_Intercept,zi_yearsSinceInf) 11588 8609
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -2.71 0.19 -3.09 -2.34 1.00 897 1550
## phi_Intercept 4.13 0.14 3.85 4.39 1.00 2833 5175
## zi_Intercept -4.77 0.73 -6.35 -3.49 1.00 3206 5944
## yearsSinceInf -0.59 0.07 -0.75 -0.45 1.00 4877 6782
## phi_yearsSinceInf 0.23 0.08 0.07 0.38 1.00 14689 9183
## zi_yearsSinceInf 0.71 0.39 -0.19 1.37 1.00 7712 6737
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# posterior predictive checks
pp_check(SuscCoral_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(SuscCoral_brm, type = "boxplot", ndraws = 6)

pp_check(SuscCoral_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(SuscCoral_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_SuscCoral*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Sponge model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_SpoCov)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_SpoCov) #-we'll work with zero-inflated beta

#define model
sponge_mod <- bf(pctCov_SpoCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = zero_inflated_beta(link = "logit", link_phi = "log", link_zi = "logit"))
##### Get priors
get_prior(sponge_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## (flat) b zi
## (flat) b yearsSinceInf zi
## logistic(0, 1) Intercept zi
## student_t(3, 0, 2.5) sd zi 0
## student_t(3, 0, 2.5) sd Site zi 0
## student_t(3, 0, 2.5) sd Intercept Site zi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site zi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_sponge <- BlackPointData %>%
filter(pctCov_SpoCov> 0)
logit(mean(BPD_sponge$pctCov_SpoCov)) #-2.4
## [1] -2.392563
##### Set priors
sponge_prior <- c(prior(normal(-2.4, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(lkj(2), class = "cor"),
prior(logistic(0, 1), class = "Intercept", dpar = "zi"),
prior(normal(0, 1), class = "b", dpar = "zi"),
prior(exponential(1), class = "sd", dpar = "zi"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"))
##### Uncomment to run model
# sponge_brm <- brm(sponge_mod,
# data = tcrmp_meta_pctcov,
# prior = sponge_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Uncomment below to save model and/or read in saved model
#saveRDS(sponge_brm, "data/outputs/sponge_brm.RDS")
sponge_brm <- readRDS("data/outputs/sponge_brm.RDS")
##### Look at model
sponge_brm
## Family: zero_inflated_beta
## Links: mu = logit; phi = log; zi = logit
## Formula: pctCov_SpoCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## zi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.70 0.09 0.54 0.91 1.00
## sd(yearsSinceInf) 0.18 0.05 0.10 0.29 1.00
## sd(phi_Intercept) 0.41 0.08 0.26 0.59 1.00
## sd(phi_yearsSinceInf) 0.28 0.11 0.07 0.52 1.00
## sd(zi_Intercept) 2.63 0.61 1.64 4.01 1.00
## sd(zi_yearsSinceInf) 0.52 0.43 0.02 1.61 1.00
## cor(Intercept,yearsSinceInf) 0.06 0.26 -0.44 0.56 1.00
## cor(phi_Intercept,phi_yearsSinceInf) 0.10 0.33 -0.51 0.75 1.00
## cor(zi_Intercept,zi_yearsSinceInf) 0.18 0.41 -0.67 0.86 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2750 4909
## sd(yearsSinceInf) 4634 7231
## sd(phi_Intercept) 5335 7016
## sd(phi_yearsSinceInf) 3521 2731
## sd(zi_Intercept) 4976 7898
## sd(zi_yearsSinceInf) 4964 6747
## cor(Intercept,yearsSinceInf) 9675 8256
## cor(phi_Intercept,phi_yearsSinceInf) 8730 7172
## cor(zi_Intercept,zi_yearsSinceInf) 15921 9099
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -3.39 0.12 -3.63 -3.15 1.00 1221 2619
## phi_Intercept 4.72 0.09 4.55 4.90 1.00 6641 7590
## zi_Intercept -5.10 0.70 -6.62 -3.90 1.00 6390 8149
## yearsSinceInf -0.04 0.05 -0.15 0.06 1.00 6727 7309
## phi_yearsSinceInf -0.02 0.11 -0.24 0.18 1.00 11167 8967
## zi_yearsSinceInf 0.43 0.41 -0.54 1.11 1.00 10812 7472
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(sponge_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(sponge_brm, type = "boxplot", ndraws = 6)

pp_check(sponge_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(sponge_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_SpoCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Turf algae model
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_EACCov)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_EACCov) #-we'll work with beta

#define model
turf_mod <- bf(pctCov_EACCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site), #-this last term allows the yearsSinceInf slope to vary across sites (also sites have varying intercepts)-#
phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf|Site),
family = Beta(link = "logit", link_phi = "log"))
##### Get priors
get_prior(turf_mod, data = tcrmp_meta_pctcov)
## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b yearsSinceInf
## lkj(1) cor
## lkj(1) cor Site
## student_t(3, 0, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd Site 0
## student_t(3, 0, 2.5) sd Intercept Site 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site 0
## (flat) b phi
## (flat) b yearsSinceInf phi
## student_t(3, 0, 2.5) Intercept phi
## student_t(3, 0, 2.5) sd phi 0
## student_t(3, 0, 2.5) sd Site phi 0
## student_t(3, 0, 2.5) sd Intercept Site phi 0
## student_t(3, 0, 2.5) sd yearsSinceInf Site phi 0
## source
## default
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## default
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
#####Calculate logit-mean of reference level for weakly-regularizing priors
BPD_turf <- BlackPointData %>%
filter(pctCov_EACCov> 0)
logit(mean(BPD_turf$pctCov_EACCov)) #-0.5
## [1] -0.5345477
##### Set priors
turf_prior <- c(prior(normal(-0.5, 1), class = "Intercept"),
prior(normal(0, 1), class = "b"),
prior(lkj(2), class = "cor"),
prior(normal(0, 1), class = "Intercept", dpar = "phi"),
prior(normal(0, 1), class = "b", dpar = "phi"),
prior(normal(0, 1), class = "sd", dpar = "phi"),
prior(exponential(1), class = "sd"))
##### Uncomment below to run model
# turf_brm <- brm(turf_mod,
# data = tcrmp_meta_pctcov,
# prior = turf_prior,
# cores = n_cores,
# chains = n_chains,
# iter = n_iter,
# warmup = n_warmup,
# init = init_list)
##### Save model and/or read in saved model
#saveRDS(turf_brm, "data/outputs/turfAlgae_brm.RDS")
turf_brm<- readRDS("data/outputs/turfAlgae_brm.RDS")
##### Look at model
turf_brm
## Family: beta
## Links: mu = logit; phi = log
## Formula: pctCov_EACCov ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## phi ~ 1 + yearsSinceInf + (1 + yearsSinceInf | Site)
## Data: tcrmp_meta_pctcov (Number of observations: 1330)
## Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 12000
##
## Group-Level Effects:
## ~Site (Number of levels: 34)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.65 0.08 0.51 0.84 1.00
## sd(yearsSinceInf) 0.16 0.04 0.09 0.25 1.00
## sd(phi_Intercept) 0.51 0.08 0.37 0.69 1.00
## sd(phi_yearsSinceInf) 0.16 0.11 0.01 0.41 1.00
## cor(Intercept,yearsSinceInf) -0.46 0.21 -0.81 -0.01 1.00
## cor(phi_Intercept,phi_yearsSinceInf) 0.02 0.38 -0.71 0.76 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2830 4796
## sd(yearsSinceInf) 4595 7342
## sd(phi_Intercept) 4712 6502
## sd(phi_yearsSinceInf) 3214 6160
## cor(Intercept,yearsSinceInf) 9183 8281
## cor(phi_Intercept,phi_yearsSinceInf) 13776 7170
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.82 0.11 -1.04 -0.60 1.00 1511 2959
## phi_Intercept 3.29 0.10 3.09 3.49 1.00 5726 7532
## yearsSinceInf 0.03 0.04 -0.06 0.11 1.00 5649 7719
## phi_yearsSinceInf 0.10 0.10 -0.09 0.31 1.00 8783 8613
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#posterior predictive checks
pp_check(turf_brm, type = "hist", ndraws = 11)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

pp_check(turf_brm, type = "boxplot", ndraws = 6)

pp_check(turf_brm, type = "dens_overlay", ndraws = 100)

y <- posterior_predict(turf_brm)*100
yrep_int <- sapply(data.frame(y), as.integer)
y_int <- as.integer(tcrmp_meta_pctcov$pctCov_EACCov*100)
ppc_bars(y_int, yrep_int)

ppc_rootogram(y_int, yrep_int)

Zoanthid model
- omitted after conversation w/ T. Smith because he said they are
extremely rare at the depths surveyed
#-Visualize changes across sites over time-#
tcrmp_meta_pctcov %>%
ggplot(aes(x = yearsSinceInf, y = pctCov_ZoanCov)) +
geom_point(aes(colour = Site)) +
geom_smooth(method = "lm", se = F, aes(colour = Site)) +
theme(legend.position = "none") +
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

#look at data distribution
hist(tcrmp_meta_pctcov$pctCov_ZoanCov) #-we'll work with zero-inflated beta

ggplot(tcrmp_meta_pctcov)+
geom_histogram(aes(x = pctCov_ZoanCov, fill = Site))+
facet_wrap(~Year)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

setup for generating predictions
#-Since our only two predictors in the model are Site and yearsSinceInf, this will let us generate predictions for all sites over a given period-#
SCTLD_sites <- unique(tcrmp_meta_pctcov$Site)
t <- seq(0, 10, by = (1/52)) #since effect is in years, this will estimate change per week
SCTLD_pred <- expand.grid(SCTLD_sites, t)
colnames(SCTLD_pred) <- c("Site", "yearsSinceInf")
Get R2 value for each model
bayes_R2(cca_brm)#.50
## Estimate Est.Error Q2.5 Q97.5
## R2 0.4992397 0.02361663 0.4518594 0.5452427
bayes_R2(cyano_brm)#.27
## Estimate Est.Error Q2.5 Q97.5
## R2 0.2652672 0.02774675 0.2129457 0.3217806
bayes_R2(gorgo_brm)#.67
## Estimate Est.Error Q2.5 Q97.5
## R2 0.6655684 0.01614981 0.6322491 0.6947094
bayes_R2(fireCoral_brm)#.24
## Estimate Est.Error Q2.5 Q97.5
## R2 0.2448967 0.03232055 0.1854114 0.3130079
bayes_R2(macro_brm)#.74
## Estimate Est.Error Q2.5 Q97.5
## R2 0.7357125 0.005802696 0.723629 0.7464166
bayes_R2(resistantCoral_brm)#.48
## Estimate Est.Error Q2.5 Q97.5
## R2 0.4765436 0.03358091 0.4115251 0.5415283
bayes_R2(SuscCoral_brm)#.80
## Estimate Est.Error Q2.5 Q97.5
## R2 0.8005005 0.007277015 0.7848155 0.8131432
bayes_R2(sponge_brm)#.62
## Estimate Est.Error Q2.5 Q97.5
## R2 0.6215961 0.01437935 0.5916809 0.6483466
bayes_R2(turf_brm)#.65
## Estimate Est.Error Q2.5 Q97.5
## R2 0.6465213 0.009391979 0.6274814 0.6641061
get posterior effect estimates for plotting
posterior_draws_cca <- cca_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "cca")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_susc <- SuscCoral_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "SCTLD-susceptible coral")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_sponge <- sponge_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "sponge")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_resistant <- resistantCoral_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "SCTLD-resistant coral")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_macroalg <- macro_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "macroalgae")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_turfalg <- turf_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "turf algae")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_gorgonian <- gorgo_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "gorgo")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_cyanobacteria <- cyano_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "cyanobacteria")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
posterior_draws_fireCoral <- fireCoral_brm %>%
posterior_samples()%>%
select(b_yearsSinceInf)%>%
rename(effect_estimate = b_yearsSinceInf)%>%
mutate(effect = "fireCoral")%>%
relocate(effect, .before = effect_estimate)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
PD_all <- posterior_draws_cca %>%
full_join(posterior_draws_cyanobacteria)%>%
full_join(posterior_draws_gorgonian)%>%
full_join(posterior_draws_macroalg)%>%
full_join(posterior_draws_resistant)%>%
full_join(posterior_draws_sponge)%>%
full_join(posterior_draws_susc)%>%
full_join(posterior_draws_turfalg)%>%
full_join(posterior_draws_fireCoral)%>%
group_by(effect)%>%
mutate(meanEffect = mean(effect_estimate))
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`
## Joining with `by = join_by(effect, effect_estimate)`